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The Lanczos algorithm for matrix tridiagonalisation suffers from strong numerical instability in finite precision 
arithmetic when applied to evaluate matrix eigenvalues. The mechanism by which this instability arises is well 
documented in the literature. A recent application of the Lanczos algorithm proposed by Bai, Fahey and Golub 
allows quadrature evaluation of inner products of the form ^ ■ g(A) ■ ip. We show that this quadrature evaluation 
is numerically stable (accurate to machine precision) and explain how the numerical errors which are such a 
fundamental element of the finite precision Lanczos tridiagonalisation procedure are automatically and exactly 
compensated in the Bai, Fahey and Golub algorithm. In the process, we shed new light on the mechanism by 
which roundoff error corrupts the Lanczos procedure. 



1. THE BAI, FAHEY AND GOLUB 
METHOD 

The determinant of the fermion matrix plays 
a central role in dynamical Lattice QCD sim- 
ulations. The direct approach jjj to evaluat- 
ing fermion determinants typically makes use of 
the Lanczos procedure for eigenvalue estimation 
of Hermitian matrices. Given an n x n Hermi- 
tian matrix A, and an orthonormal seed vector 
/i, the Lanczos procedure is an iterative process 
which (in exact arithmetic) generates a sequence 
of orthogonal vectors fi, i = 1, . . . , m. If we de- 
fine an n x m matrix F m = (/i, . . . , f m ) whose 
columns are given by the vectors /,;, then we have: 
F L F m = l m , T m = F+jAF m is tridiagonal, and the 
eigenvalues of T m converge to the eigenvalues of 
A as m — > n. 

An alternative stochastic approach to evaluat- 
ing determinants was proposed in 1996 by Bai, 
Fahey and Golub (BFG) P] who observe that a 
vector-matrix- vector of the form ip'-g(A)-ip can be 
expressed as an integral of the function g(-) over a 
modified spectral measure. An m-point Gaussian 
quadrature approximation to this integral is then 
given by xjfi -g(A)-ip sa Y^Li w i9(@i)- The remark- 
able result J! applied in the BFG method is that 



the abscissa, Qi , of this quadrature rule are given 
by the eigenvalues of the tridiagonal matrix T m 
which is generated by a Lanczos procedure ap- 
plied with seed vector ip, and that the weights, 
Wi, are given as the squares of the first compo- 
nents of the corresponding eigenvectors of T m . 

Figure [l] shows an example of the results which 
the BFG method generates on a 12 3 x 24 quenched 
gauge configuration. The configuration is a 
quenched configuration with f3 = 5.7 produced 
within the current UKQCD research program. A 
random Gaussian fermion vector, ip, was gener- 
ated, and the BFG method applied to evaluate 
ip' ■ log(Mj_M K ) • ip for Wilson hopping parame- 
ter k = .1650 as a function of the order m of the 
Lanczos tridiagonal matrix and of the resulting 
Gaussian quadrature approximation. To generate 
an estimate of the "true" value for this matrix- 
vector inner-product for comparison, we arbitrar- 
ily averaged the values for order 1991 to order 
2000. The plot then shows the log of the differ- 
ence of the m-th order BFG estimate and this 
"true" value plotted against m. The plot shows 
a rapid convergence to the "true" value (conver- 
gence has occurred at about order 200), followed 
by what is clearly fluctuations at the level of algo- 
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Figure 1. Convergence history for the evaluation 
of ip^ -log(Mj.M K ) on a 12 3 x 24 quenched con- 
figuration as a function of the number of Lanczos 
iterations. 



rithm precision. Note in particular that the long 
converged region is perfectly flat, and shows no 
drift. Checks for a toy model, and the mathe- 
matical analysis which we briefly describe here, 
indicate that the convergence is in fact, conver- 
gence to the exact answer, and that finite preci- 
sion arithmetic introduces no systematic error. 

2. THE EFFECT OF FINITE PRECI- 
SION ARITHMETIC 

At first sight, Figure [I] might be interpreted as 
good phenomenological evidence that the BFG 
method is reliable and robust. A second sight 
however should disturb a person familiar with 
the Lanczos procedure at finite precision. In fi- 
nite precision arithmetic, the Lanczos procedure 
generates clusters of almost degenerate eigenval- 
ues of the tridiagonal matrix T m where only a 
single eigenvalue would have occurred in an in- 
finite precision calculation. These clusters could 
potentially destroy the convergence of the BFG 
quadrature sum to the correct infinite precision 
result since the corresponding abscissa are in- 
cluded multiple times rather than singly in that 
sum. We have found however, that when clusters 
arise, the corresponding weights in the quadra- 



ture rule adjust so that the net contribution of 
the cluster of multiple eigenvalues and weights is 
equal to that of the corresponding single eigen- 
value and weight which would arise in infinite 
precision. 

This effect is illustrated in Table |l| which shows 
the weights corresponding to a given cluster for 
a Lanczos procedure applied to a model problem 
which has been chosen so all quantities can be cal- 
culated both numerically and exactly and direct 
comparisons can be made. At iteration 40, the 
cluster in question contains only a single eigen- 
value, and the relative error between the value of 
the weight as generated by the Lanczos procedure 
and an exact calculation is seen to be O(10~ 14 ). 
By iteration 160, the cluster expands to include 
four almost degenerate eigenvalues (which typi- 
cally differ only at about O(10~ 8 )). The individ- 
ual weights are seen now to differ quite signifi- 
cantly from the correct value, but the sum of the 
three weights does match the correct value, and 
the contribution of the cluster to the quadrature 
is seen to be independent of the size of the cluster. 

3. AN EXPLANATION 

The eigenvalue clustering effect which is well 
known J4|, and the invariance of weights which 
we have just described imposes strict constraints 
on the mechanism by which roundoff error cor- 
rupts the Lanczos procedure. At finite precision, 
the Lanczos procedure still generates a sequence 
of vectors fa, i = 1 . . . m. However these vectors 
no longer remain orthogonal, and no longer tridi- 
agonalise A. We have instead that 

2. Fj^AFm — H m , symmetric but not tridiago- 
nal, 

3. The tridiagonal matrix, T m which is, in 
practice, used to estimate eigenvalues of A is 
given by the truncation of H m to tridiagonal 
form (i.e. the all off-tridiagonal elements of 
H m are set to zero). 

Ignoring the very small differences be- 
tween eigenvalues in a degenerate cluster, 
the eigenvalue spectrum of T m is given as 
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Table 1 

This table illustrates the invariance of the sum of the weights for an eigenvalue cluster generated from a 
finite precision Lanczos procedure. 



m Weights in Cluster 



40 
120 



160 



Sum of Weights 



Relative Error 



0.121274864637753313E-02 
0.239992265988143331E-04 
0.298384753074422057E-03 
0.890364666704296305E-03 
0.128211917795651532E-05 
0.1 18900026992380976E-03 
0.470001178219875870E-03 
0.622565321987319131E-03 



0.121274864637753313E-02 0.250E-14 
0.121274864637753270E-02 0.286E-14 



0.121274864637753248E-02 0.303E-14 



{6i(ni),...,9j(nj),9j+i(l),...9 p (l)}, where m 
is the multiplicity of the cluster with eigenvalue 
0i, j counts the non- unit-multiplicity clusters, 
and p counts the different distinct eigenvalues. 
We denote the rij eigenvectors with degenerate 

(k) 

eigenvalue 9i by si ,k = 1 . . . n% . The n-row 
by m column matrix F m relates these eigenvec- 
tors to "Ritz" vectors yi which are n-component 
vectors in the space spanned by the columns of 

(k) 

F m as yi = F m s- . A well established result 

H for the Lanczos procedure is that, in a non- 

(k) 

unit-multiplicity cluster, the eigenvectors s] all 
generate the same Ritz vector yi (approximately). 
Further this Ritz vector is a true eigenvector of 
the matrix A, and the corresponding eigenvalue 
is a true eigenvalue of A (approximately). The 
matrix F m is therefore seen to be rank-deficient 
since there are independent linear combinations 
of columns which equal the same vector. 

We have found || that the best way to ex- 
pose this rank deficiency is to execute a singu- 
lar value decomposition of F m . Using this singu- 
lar value decomposition, it is possible to find a 
special set of orthonormal basis vectors for the 
space spanned by the columns of F m . In this 
special basis, the matrix H m takes a block diag- 
onal form with blocks corresponding to the non- 
unit-multiplicity eigenvalues of T m having entries 
at all positions equal to the corresponding eigen- 
value. When H m is truncated to diagonal form in 
the normal way, these blocks in the special basis 
representation become diagonal (the effect of the 
truncation is to set all the non-diagonal entries for 



the non- unit-multiplicity blocks equal to zero). 
We thus obtain the usual degenerate eigenvalue 
clusters of the truncated matrix T m . The reason 
that the Lanczos procedure at finite precision pro- 
duces good eigenvalue estimates is thus seen to be 
a combination of the specially structured nature 
of the space spanned by the columns of F m , and 
of the very special way in which the truncation to 
diagonal form acts. The invariance of the sum of 
weights for a cluster is then explained by the ob- 
servation that this sum gives the projection of the 
starting vector ip seeding the Lanczos procedure 
onto the corresponding Ritz vector. 
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